话说当年我在暨大工作的时候，学生都很害怕汇报RT-PCR的结果，因为每次一讲到用了t test做了检验，底下某小老板，就一顿狂喷，说学生瞎搞，这不对那不对，但从没告诉正确的姿势，学生们都怕了他。于是我写了这篇文章，普及了一下相关的知识，RT-PCR那点东西，全在这里。
RT-PCR相对定量有两个模型，一是ΔΔCt 模型，一个是扩增效率校正模型。这里我们先讨论简单的模型：ΔΔCt 模型，在这一模型中，假定扩增效率为2，即每个PCR cycle，产物倍增，由以下公式给出：
Ratio=2−ΔΔCt
其中ΔΔCt = ΔΔCtreated - ΔΔCcontrol, ΔΔCtreated和ΔΔCcontrol分别是treatment组和control组中目标基因和参照基因的Ct差，即ΔΔCt = ΔΔCtarget - ΔΔCreference.
扩增效率肯定是在(1,2)之间，理想状态下才是2，而理想状态通常是不存在的。这个模型的假设点其实是treatment组和control组的扩增效率一样，至于是不是2，比2小多少，并不影响后续的统计分析。
Data
我们来看以下一份数据，4组重复实验，用RT-PCR测了gene01和gene02的表达量，HK代表house keeping gene，即参照。以下数值为原始的Ct值。
ct <- data.frame(
sample = rep(rep(1:4, each=3), 2),
treatment = rep(c("Control","Treated"), each=12),
gene01 = c(23.22,23.34,23.12,24.06,24.15,24.15,23.18,23.13,
23.10,24.78,24.45,24.67,23.11,22.99,23.10,22.77,
22.99,23.06,23.73,24.01,23.80,23.73,23.83,23.73),
gene02 = c(29.08,29.04,29.39,28.23,28.01,28.12,28.79,28.43,
28.49,31.37,30.74,31.09,27.11,27.24,27.37,25.52,
25.72,25.52,27.43,26.73,26.65,27.96,27.84,27.98),
HK = c(19.68,19.69,19.80,19.95,19.93,19.97,19.93,20.02,20.27,
19.93,19.88,19.90,20.61,19.98,20.57,19.68,19.95,
19.85,20.27,20.08,20.07,20.10,20.07,20.10)
)
print(ct)
# sample treatment gene01 gene02 HK
# 1 1 Control 23.22 29.08 19.68
# 2 1 Control 23.34 29.04 19.69
# 3 1 Control 23.12 29.39 19.80
# 4 2 Control 24.06 28.23 19.95
# 5 2 Control 24.15 28.01 19.93
# 6 2 Control 24.15 28.12 19.97
# 7 3 Control 23.18 28.79 19.93
# 8 3 Control 23.13 28.43 20.02
# 9 3 Control 23.10 28.49 20.27
# 10 4 Control 24.78 31.37 19.93
# 11 4 Control 24.45 30.74 19.88
# 12 4 Control 24.67 31.09 19.90
# 13 1 Treated 23.11 27.11 20.61
# 14 1 Treated 22.99 27.24 19.98
# 15 1 Treated 23.10 27.37 20.57
# 16 2 Treated 22.77 25.52 19.68
# 17 2 Treated 22.99 25.72 19.95
# 18 2 Treated 23.06 25.52 19.85
# 19 3 Treated 23.73 27.43 20.27
# 20 3 Treated 24.01 26.73 20.08
# 21 3 Treated 23.80 26.65 20.07
# 22 4 Treated 23.73 27.96 20.10
# 23 4 Treated 23.83 27.84 20.07
# 24 4 Treated 23.73 27.98 20.10
Data clean
每个sample，测了三次技术重复，我们使用平均值来提高精度。
require(plyr)
# Loading required package: plyr
ct <- ddply(ct, .(sample, treatment), function(x)
data.frame(gene01=mean(x$gene01),
gene02=mean(x$gene02), HK=mean(x$HK)))
print(ct)
# sample treatment gene01 gene02 HK
# 1 1 Control 23.22667 29.17000 19.72333
# 2 1 Treated 23.06667 27.24000 20.38667
# 3 2 Control 24.12000 28.12000 19.95000
# 4 2 Treated 22.94000 25.58667 19.82667
# 5 3 Control 23.13667 28.57000 20.07333
# 6 3 Treated 23.84667 26.93667 20.14000
# 7 4 Control 24.63333 31.06667 19.90333
# 8 4 Treated 23.76333 27.92667 20.09000
Calculate statistical metric
按照ΔΔCt 模型，我们先计算ΔCt，然后计算ΔΔCt，最终计算出Ratio。
delta.ct <- ddply(ct, .(sample, treatment), function(x)
data.frame(gene01=x$gene01-x$HK,
gene02=x$gene02-x$HK))
print(delta.ct)
# sample treatment gene01 gene02
# 1 1 Control 3.503333 9.446667
# 2 1 Treated 2.680000 6.853333
# 3 2 Control 4.170000 8.170000
# 4 2 Treated 3.113333 5.760000
# 5 3 Control 3.063333 8.496667
# 6 3 Treated 3.706667 6.796667
# 7 4 Control 4.730000 11.163333
# 8 4 Treated 3.673333 7.836667
# Delta Delta Ct
dd.ct <- subset(delta.ct, treatment=="Treated",
select=c(gene01, gene02)) -
subset(delta.ct,
treatment=="Control",
select=c(gene01, gene02))
print(dd.ct)
# gene01 gene02
# 2 -0.8233333 -2.593333
# 4 -1.0566667 -2.410000
# 6 0.6433333 -1.700000
# 8 -1.0566667 -3.326667
ratio <- 2^(-1*dd.ct)
print(ratio)
# gene01 gene02
# 2 1.769490 6.034915
# 4 2.080120 5.314743
# 6 0.640232 3.249010
# 8 2.080120 10.032899
这个ratio值计算出，经过扩增后，目的基因在处理组中的产量是对照组的多少倍。
Statistical Analysis
计算出ratio后，很多人就拿它来做t检验。这种做法是不对的，因为ratio的分布明显是右偏的。从理论上看，PCR的扩增可以分成以下三个阶段，指数扩增、线性扩增、平台期。
当试剂量足够的时候，就处于指数扩增阶段，只有当PCR试剂不足时，才进入线性扩增期，最后当试剂被耗尽时，产物不增长，位于平台期。
在我们的实验中，通常PCR反应是处于指数增长期。我们将其进行对数处理，指数增长就会变成线性增长，这样数据的右偏现象就会有明显减弱，当然我不敢说对数转换后，它就能呈正态分布，但起码对数处理后数据的分布比原始ratio数据更接近正态了。
T Test
只有偏态不明显的情况下，才可以用t检验进行分析。
# T test for gene01
t.test(log(ratio[,1], base=2), mu = 0)
#
# One Sample t-test
#
# data: log(ratio[, 1], base = 2)
# t = 1.4009, df = 3, p-value = 0.2558
# alternative hypothesis: true mean is not equal to 0
# 95 percent confidence interval:
# -0.729139 1.875806
# sample estimates:
# mean of x
# 0.5733333
# T test for gene02
t.test(log(ratio[,2], base=2), mu = 0)
#
# One Sample t-test
#
# data: log(ratio[, 2], base = 2)
# t = 7.5039, df = 3, p-value = 0.004904
# alternative hypothesis: true mean is not equal to 0
# 95 percent confidence interval:
# 1.44405 3.57095
# sample estimates:
# mean of x
# 2.5075
Use −ΔΔCt directly
如果我们要使用T检验，就得做对数处理，细心的人应该会发现ratio的计算就是指数运算，ratio取对数之后，就是−ΔΔCt，所以我们可以直接使用它来做检验。
# T test for gene02
t.test(-dd.ct[,2])
#
# One Sample t-test
#
# data: -dd.ct[, 2]
# t = 7.5039, df = 3, p-value = 0.004904
# alternative hypothesis: true mean is not equal to 0
# 95 percent confidence interval:
# 1.44405 3.57095
# sample estimates:
# mean of x
# 2.5075
ANOVA
另一种方法，使用ANOVA进行统计，这也是参数统计方法，在我们这种小样本的情况下，同样对数据分布有严格的要求，我们不能基于ratio做统计，要基于ΔCt值。
print(delta.ct)
# sample treatment gene01 gene02
# 1 1 Control 3.503333 9.446667
# 2 1 Treated 2.680000 6.853333
# 3 2 Control 4.170000 8.170000
# 4 2 Treated 3.113333 5.760000
# 5 3 Control 3.063333 8.496667
# 6 3 Treated 3.706667 6.796667
# 7 4 Control 4.730000 11.163333
# 8 4 Treated 3.673333 7.836667
summary(aov(gene01 ~ treatment, data=delta.ct))
# Df Sum Sq Mean Sq F value Pr(>F)
# treatment 1 0.6574 0.6574 1.687 0.242
# Residuals 6 2.3385 0.3898
summary(aov(gene02 ~ treatment, data=delta.ct))
# Df Sum Sq Mean Sq F value Pr(>F)
# treatment 1 12.575 12.575 9.963 0.0197 *
# Residuals 6 7.573 1.262
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ANOVA也是比较常用的，而且可以应用于多组实验和多个因素中。但是如果按照上面的代码来做，还不如用T test，因为实验是成对（tretated vs control）进行的，而这个信息在上面的ANOVA分析中，被扔了。
实际上 −ΔΔCt 的单样本T检验，相当于ΔΔCt 的双样本T检验。
t.test(gene02 ~ treatment, delta.ct, paired=T)
#
# Paired t-test
#
# data: gene02 by treatment
# t = 7.5039, df = 3, p-value = 0.004904
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# 1.44405 3.57095
# sample estimates:
# mean of the differences
# 2.5075
而这里用ANOVA，没有用到成对的信息，计算出来的p值和unpaired T test就基本相等。
t.test(gene01 ~ treatment, delta.ct, paired=F)
#
# Welch Two Sample t-test
#
# data: gene01 by treatment
# t = 1.2988, df = 5.2396, p-value = 0.2482
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# -0.5460188 1.6926855
# sample estimates:
# mean in group Control mean in group Treated
# 3.866667 3.293333
t.test(gene02 ~ treatment, delta.ct, paired=F)
#
# Welch Two Sample t-test
#
# data: gene02 by treatment
# t = 3.1565, df = 5.064, p-value = 0.02476
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# 0.473167 4.541833
# sample estimates:
# mean in group Control mean in group Treated
# 9.319167 6.811667
做为成对T检验的扩展，使用方差分析进行多组实验的比较，需要用到的是重复测量方差分析（repeated-measures ANOVA）。
delta.ct[,1] = factor(delta.ct[,1])
delta.ct[,2] = factor(delta.ct[,2])
summary(aov(gene01 ~ treatment+Error(sample/treatment), data=delta.ct))
#
# Error: sample
# Df Sum Sq Mean Sq F value Pr(>F)
# Residuals 3 1.333 0.4445
#
# Error: sample:treatment
# Df Sum Sq Mean Sq F value Pr(>F)
# treatment 1 0.6574 0.6574 1.962 0.256
# Residuals 3 1.0050 0.3350
summary(aov(gene02 ~ treatment+Error(sample/treatment), data=delta.ct))
#
# Error: sample
# Df Sum Sq Mean Sq F value Pr(>F)
# Residuals 3 6.903 2.301
#
# Error: sample:treatment
# Df Sum Sq Mean Sq F value Pr(>F)
# treatment 1 12.57 12.575 56.31 0.0049 **
# Residuals 3 0.67 0.223
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
这样算出来的p值就和成对T检验几乎一样。
Wilcoxon Rank Sum Test
我们计算出来的ratio，是经过PCR扩增之后产物的比值，这是经过指数扩增的，而如果我们取对数，相当于消除了指数增长的效果，−ΔΔCt 其实就是原始量的差值。
所以ratio的比较，相当于是扩增后的比较，而log2ratio 的比较，则相当于是扩增之前的比较。这两种数据当然都是可以拿来做统计分析的，只不过T test等参数检验方法对数据分布有要求，不适合于ratio值的统计分析而已。如果我们使用非参数检验，则不管是ratio值也好，log2ratio值也好，结果都是一样的。非参数检验并不使用数据的原始值，而是使用其排序值来做检验，我们并不能确定对数处理后的数据就呈正态分布，还有就是通常我们的实验重复次数是非常少的，在这种情况下，使用非参数检验，其实也是一种保守的、不容易犯错的方法。
# Wilcoxon Rank Sum test for gene02, original data
wilcox.test( ratio[,2], mu=1)
#
# Wilcoxon signed rank test
#
# data: ratio[, 2]
# V = 10, p-value = 0.125
# alternative hypothesis: true location is not equal to 1
# Wilcoxon Rank Sum test for gene02, log2 transformed data
wilcox.test( log(ratio[,2], base=2), mu=0)
#
# Wilcoxon signed rank test
#
# data: log(ratio[, 2], base = 2)
# V = 10, p-value = 0.125
# alternative hypothesis: true location is not equal to 0
说它不容易犯错，是因为它对数据分布没要求；说它保守，是因为在计算p-value上，灵敏度没有T test这么高，p-value会大一些。
Permutation Test
生物学的数据通常比较messy，non-normal，重复次数又少。在我们不能确定数据分布是否正态的情况下，使用参数检验，很可能结果bias太大，使用非参数检验，又觉得power不够大，那么我们可以使用permutation test。这不失为一种好方法，通过随机变换两组数据的标签，观察它们的均值差，我们将得到均值差的分布，将我们实际的测量值对比于随机得到的差值分布，算出随机得到的差值大于实际的测量值的概率。
perm.test <- function(d1, d2, nIter=10000) {
if (length(d2) == 1)
d2 <- rep(d2, length(d1))
m <- mean(d2-d1)
pooledData <- c(d1, d2)
n <- length(d1)
meanDiff <- numeric(nIter)
for (i in 1:nIter) {
idx <- sample(1:length(pooledData), n, replace=FALSE)
d1 <- pooledData[idx]
d2 <- pooledData[-idx]
meanDiff[i] <- mean(d2) - mean(d1)
}
p <- mean(abs(meanDiff) >= abs(m))
return(p)
}
perm.test(ratio[,1],1)
## [1] 0.1388
perm.test(-dd.ct[,1], 0)
## [1] 0.1454
perm.test(ratio[,2],1)
## [1] 0.0288
perm.test(-dd.ct[,2], 0)
## [1] 0.0266
从结果上看，用ratio还是用−ΔΔCt 差不多，这表明原始数据的分布完全不关事。
没有吃没有穿
自有那敌人送上前
没有枪没有炮
敌人给我们造
----游击队之歌
这个方法的奇妙之处，就在于通过bootstrap，产生针对我们统计量的一个背景分布，再由这个背景分布计算出p-value。
Plot
至于画图，我推荐使用−ΔΔCt 来画，这个CONTROL是0，不用画。如果要用ratio来画，你觉得CONTROL每次都是画一个mean=1，sd=0的柱子，有意义吗？anyway，我下面的画图代码，同样也适用于拿来画ratio值。
BTW：这里图中标的p值，使用的是T检验算出来的p值。
require(reshape2)
## Loading required package: reshape2
dd.ctm <- melt(-dd.ct)
## No id variables; using all as measure variables
print(dd.ctm)
## variable value
## 1 gene01 0.8233333
## 2 gene01 1.0566667
## 3 gene01 -0.6433333
## 4 gene01 1.0566667
## 5 gene02 2.5933333
## 6 gene02 2.4100000
## 7 gene02 1.7000000
## 8 gene02 3.3266667
dd.cts <- ddply(dd.ctm, .(variable), function(x)
data.frame(mean=mean(x$value), sd=sd(x$value)))
print(dd.cts)
## variable mean sd
## 1 gene01 0.5733333 0.8185353
## 2 gene02 2.5075000 0.6683222
require(ggplot2)
## Loading required package: ggplot2
p <- ggplot(dd.cts, aes(variable, mean, fill=variable, width=.5))+
geom_bar(stat="identity", position=position_dodge(width=.8), colour="black")+
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2, position=position_dodge(width=.8))
p <- p+scale_fill_manual(values=c("white", "grey80"))
p <- p+theme_bw()+theme(legend.position="none") + xlab("")+ylab(expression(paste("-", Delta, Delta, "Ct")))
p <- p+theme(axis.text.x=element_text(face="bold", size=14), axis.text.y=element_text(face="bold", size=14), axis.title.y=element_text(size=18, face="bold"))
p <- p+annotate("text", x=2,y=3.5, label="p < 0.005\n**", size=4)
print(p)
再来一个箱式图，我其实觉得箱式图比柱状图要好一些，这里我除了用箱式图画数据分布之后，把原始的数据点也给画了。
p2 <- ggplot(dd.ctm, aes(variable, value, fill=variable, colour=variable))
p2 <- p2+geom_boxplot(alpha=.3, width=.5)+geom_dotplot(binaxis='y', stackdir='center', dotsize=.5)
p2 <- p2+theme_bw()+theme(legend.position="none") + xlab("")+ylab(expression(paste("-", Delta, Delta, "Ct")))
p2 <- p2+theme(axis.text.x=element_text(face="bold", size=14), axis.text.y=element_text(face="bold", size=14), axis.title.y=element_text(size=18, face="bold"))
p2 <- p2+annotate("text", x=2,y=3.5, label="p < 0.005\n**", size=4)
print(p2)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
Summary
据我的观察，生物坑里大家对T检验特别有爱，我甚至看到用T检验去分析免疫组化数据的，而且还发在Nature子刊上，我都吐槽无力了。如果你喜欢T检验，你要注意数据分布，那么在分析RT-PCR数据上，请用−ΔΔCt作为统计量，或者将你算完的ratio值取对数也是一样的。用ANOVA时，也必须同样注意。
在这个重复次数少，数据又非常messy的生物世界里，保守点的做法是使用非参数统计，即使你使用ratio这种严重右偏的数据，你的统计分析照样没问题，多学点总是好的，凡事不能太依赖于T检验，因为在很多其它的场景里，正如这里的ratio值，你不能简单套用T检验。
Permutation test其实才是工具箱里不可或缺的一件，应用场景和非参一样广，同时灵敏度又比非参高，你值得拥有。